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We develop numerical schemes for solving the isothermal compressible and incompressible 
equations of fluctuating hydrodynamics on a grid with staggered momenta. We develop a 
second-order accurate spatial discretization of the diffusive, advective and stochastic fluxes 
that satisfies a discrete fluctuation-dissipation balance, and construct temporal discretiza- 
tions that are at least second-order accurate in time deterministically and in a weak sense. 
Specifically, the methods reproduce the correct equilibrium covariances of the fiuctuating 
fields to third (compressible) and second (incompressible) order in the time step, as we ver- 
ify numerically. We apply our techniques to model recent experimental measurements of 
giant fiuctuations in diffusively mixing fluids in a micro-gravity environment [A. Vailati et. 
ai. Nature Communications 2:290, 2011]. Numerical results for the static spectrum of non- 
equilibrium concentration fluctuations are in excellent agreement between the compressible 
and incompressible simulations, and in good agreement with experimental results for all 
measured wavenumbers. 

I. INTRODUCTION 

At a molecular scale, fluids are not deterministic; the state of the fluid is constantly changing 
and stochastic, even at thermodynamic equilibrium. Stochastic effects are important for flows 
in new microfluidic, nanofluidic and microelectromechanical devices jl]; novel materials such as 
nanofluids [2]; biological systems such as lipid membranes [3], Brownian molecular motors [1], 
nanopores [5j; as well as processes where the effect of fluctuations is amplified by strong non- 
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equilibrium effects, such as combustion of lean flames, capillary dynamics [6", T], and hydrodynamic 
instabilities [HHlOj. and others. Because they span the whole range of scales from the microscopic 
to the macroscopic |1H [T^ , fluctuations need to be consistently included in all levels of description 
|13j . Thermal fluctuations are included in the fluctuating Navier-Stokes (NS) equations and related 
continuum Langevin models [141 115j through stochastic forcing terms, as first proposed by Landau 
and Lifshitz [161. Numerically solving the continuum equations of fluctuating hydrodynamics |17| 
is difficult because of the presence of non-trivial dynamics at all scales and the existence of a 
nontrivial invariant measure (equilibrium distribution). 

Several numerical approaches for fluctuating hydrodynamics have been proposed. The earli- 
est work by Garcia et al. |18j developed a simple scheme for the stochastic heat equation and 
the linearized one-dimensional fluctuating NS equations. Ladd and others have included stress 
fluctuations in (isothermal) Lattice Boltzmann methods for some time [19j. Moseler and Land- 
man [8J included the stochastic stress tensor of Landau and Lifshitz in the lubrication equations 
and obtain good agreement with their molecular dynamics simulation in modeling the breakup of 
nanojets. Sharma and Patankar [20j developed a fluid-structure coupling between a fluctuating in- 
compressible solver and suspended Brownian particles. Coveney, De Fabritiis, Delgado-Buscalioni 
and co-workers have also used the fluctuating isothermal NS equations in a hybrid scheme, coupling 
a continuum fluctuating solver to a molecular dynamics simulation of a liquid |21H23j . Atzberger, 
Kramer and Peskin have developed a version of the immersed boundary method that includes 
fluctuations j241 125j. Voulgarakis and Chu [26] developed a staggered scheme for the isothermal 
compressible equations as part of a multiscale method for biological applications, and a similar 
staggered scheme was also described in Ref. [27j. 

Some of us have recently developed techniques for analyzing the weak accuracy of finite- volume 
methods for solving the types of stochastic partial differential equations that appear in fluctuat- 
ing hydrodynamics [2H]. The analysis emphasizes the necessity to maintain fluctuation-dissipation 
balance in spatio-temporal discretizations [28] , thus reproducing the Gibbs-Boltzmann distribution 
dictated by equilibrium statistical mechanics. Based on previous work by Bell et al. |29l I30j . a 
collocated spatial discretization for the compressible equations of fluctuating hydrodynamics has 
been developed and combined with a stochastic third-order Runge-Kutta (RK3) temporal integra- 
tor [28]. The collocated spatial discretization has been used to construct a strictly conservative 
particle-continuum hybrid method [TH] and to study the contribution of advection by thermal 
velocities to diffusive transport |31| . 

A staggered spatial discretization is advantageous for incompressible flows because it leads to 
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a robust idempotent discrete projection operator [32l [33]. Staggered schemes have previously 
been developed for isothermal compressible [26] and incompressible flow |20j . without, however, 
carefully assessing discrete fluctuation-dissipation balance. Here we present and test an explicit 
compressible and a semi-implicit incompressible scheme for fluctuating hydrodynamics on uniform 
staggered grids. Both methods use closely-related spatial discretizations, but very different tempo- 
ral discretizations. In the spatial discretization, we ensure an accurate spectrum of the steady-state 
fluctuations by combining a locally-conservative finite- volume formulation, a non-dissipative (skew- 
symmetric) advection discretization, and discretely dual divergence and gradient operators. For 
compressible flow, we employ an explicit RK3 scheme |28) since the time step is limited by the 
speed of sound and the dissipative terms can be treated explicitly. For incompressible flow, we 
use a semi-implicit unsplit method first proposed in Ref. [34j, which allows us to take large time 
steps that under-resolve the fast momentum diffusion at grid scales but still obtain the correct 
steady-state covariances of fiuctuations. 




Figure 1: Snapshots of concentration showing the development of a rough diffusive interface between two 
miscible fluids in zero gravity. We show three points in time (top to bottom), starting from an initially 
perfectly flat interface (phase separated system) . These figures were obtained using the incompressible code 
described in Section HV Al 

Thermal fiuctuations in non-equilibrium systems in which a constant (temperature, concentra- 
tion, velocity) gradient is imposed externally exhibit remarkable behavior compared to equilibrium 
systems. Most notably, external gradients can lead to enhancement of thermal fluctuations and 
to long-range correlations between fluctuations |17[ ESHSH] • This phenomenon can be illustrated 
by considering concentration fluctuations in an isothermal mixture of two miscible fluids in the 
presence of a strong concentration gradient Vc, as in the early stages of diffusive mixing between 
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initially separated fluid components. As illustrated in Fig. [T| the interface between the fluids, 
instead of remaining flat, develops large-scale roughness that reaches a pronounced maximum until 
gravity or boundary effects intervene. These giant fluctuations |39I - HT] during free diffusive mixing 
have been observed using light scattering and shadowgraphy techniques \12\ H2l - I45j , finding good 
but imperfect agreement between the predictions of a simplified fluctuating hydrodynamic theory 
and experiments. Recent experiments have taken advantage of the enhancement of the nonequi- 
librium fluctuations in a microgravity environment aboard the FOTON M3 spaceship |12l 144] , and 
demonstrated the appearance of fractal diffusive fronts like those illustrated in Fig. [1} In the ab- 
sence of gravity, the density mismatch between the two fluids does not change the qualitative nature 
of the non-equilibrium fluctuations, and in this work we focus on mixtures of dynamically-identical 
fluids. 

Before discussing spatio-temporal discretizations, we review the continuum formulation of the 
equations of fluctuating hydrodynamics and their crucial properties in Section In particular, we 
discuss the steady-state covariances of the fluctuating fields for systems in thermal equilibrium as 



well as fluid mixtures with an imposed concentration gradient. In Section III A we focus on the 
temporal discretization in the spirit of the method of lines. For the compressible equations, we 
employ a previously-developed explicit three-stage Runge-Kutta scheme that is third order weakly 
accurate p8]. For the incompressible equations we employ a second-order accurate predictor- 
corrector approach, each stage of which is a semi-implicit (Crank-Nicolson) discretization of the 
Stokes equations, solved effectively using a projection method as a preconditioner |34j . In Section 



IIIB5 we describe a conservative staggered spatial discretization of the diffusive, stochastic and 



advective fluxes. We maintain discrete fluctuation-dissipation balance |28[ 146] by ensuring duality 
between the discrete divergence and gradient operators, and by using a skew-adjoint discretization 
of advection. We verify the weak order of accuracy for both the compressible and incompressible 



algorithms in Section IV In Section |V] we model the non-equilibrium concentration fluctuations 
in a fluid mixture under an applied temperature gradient, and compare the numerical results to 
recent experimental measurements [121 



A. Deterministic Hydrodynamic Equations 

Motivated by the microgavity experiments studying giant fluctuations |12l I44j . we consider an 
ideal solution of a macromolecule with molecular mass M in a solvent. At the macroscopic level, the 
hydrodynamics of such a mixture can be modeled with an extended set of Navier-Stokes equations 
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for the mass density p = pi+ P2, where pi is the mass density of the solute, v is the center-of-mass 
velocity, c = pi/p is the mass concentration, and T is the temperature |17l I30j . In many situations 
of interest the temperature T {r,t) = T (r) can be taken as fixed \n\ \W\ BS] since temperature 
fluctuations do not couple to other variables. The fixed-temperature compressible NS equations 
for an ideal mixture of two miscible fluids are 

Dtp = -piV-v) (1) 
p{Dtv) = - VP + V • [rfVv + C (V • /] + (2) 
p {Dtc) = V • [px (Vc + c (1 - c) SrVT)] + /„ (3) 

supplemented with appropriate boundary conditions. Here DfD = SjD + • V (□) is the advective 
derivative, = (Vv + Vt?"^) — 2 ( V • v)I/3 is the symmetrized strain rate, P{p,c;T) is the 
pressure as given by the equation of state, and and fc are external forcing (source) terms. The 
shear viscosity r/, bulk viscosity C,, mass diffusion coefficient X: and Soret coefficient St, can, in 
general, depend on the state. We make several physically-motivated approximations, including 
neglecting barodiffusion, as we describe and justify next. 

We will assume that the two species in the mixture are almost identical, meaning that none of 
the fluid properties are affected by concentration. In this sense the macromolecules are assumed to 
be simple (passive) tracer particles. This is a reasonable approximation for small concentrations 
c ^ 1, since the presence of small amounts of macromolecule causes small changes in the properties 
of the solution. In the giant fluctuation experiments in microgravity conditions the concentration is 
at most a few percent \12\ I44j. justifying the assumption that the equation of state is independent 
of concentration, P{p,c;T) = P{p;T). This approximation also allows us to neglect barodiffu- 
sion since the barodiffusion coefficient is a thermodynamic rather than a transport coefficient and 
vanishes for such an equation of state. 

Because the temperature varies by only a few percent across the sample in the giant fluctuation 
experiments modeled in Section |Vj we take the system to be isothermal and thus the temperature 
T = To is constant. However, we retain the crucial Soret term by taking S'^VT to be a specified 
constant. Note that the Soret term is a transport coefficient unlike the barodiffusion coefficient 
and can be positive or negative. 

For liquids, the equation of state is usually very stiff, which means that the (isothermal) sound 
speed = dP/dp is very large. Density is therefore nearly constant, and an incompressible 
approximation will be appropriate as a means to avoid the stiffness. An alternative, employed 
for example in the Lattice-Boltzmann method, is to keep the simpler compressible equations (and 
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thus avoid elliptic constraints), but make the speed of sound much smaller than the actual speed 
of sound, but still large enough that density variations are negligible. This is the sense in which 
we will use the compressible equations 
which it is actually important to solve these equations with the proper equation of state [261 147j. 

Under the assumption that density variations are small, it is not important what precise depen- 
dence of the transport coefficients and equation of state on the density is used. We will therefore 
assume that P = P{p) = Pq + [p — po) c^, where ct is a spatially-constant isothermal speed of 
sound. The value of ct can be a parameter that lets us tune the compressibility or the physical 
speed of sound. Furthermore, we will assume that the viscosity and Soret coefficient are constants 
independent of the density, and that the product px = PoXo is constant. Recalling that in the 
experiments c <C 1 so that c {1 — c) ^ c, all of these approximations allows us to write the viscous 
term in the momentum equation in the "Laplacian" form 

V • [7]Vv + C{V -v)!]^ r]V^v + (c + ^) V (V • . (4) 

Similarly, the diffusive term in the concentration equation can be written as 

V • [px (Vc + c (1 - c) StVT)] ^ p [xV^c + V • (c^,)] , (5) 

where the spatially-constant velocity difference between the two species is denoted with Vg = 
xStVT. 

With all of these simplifications, the equations we actually solve numerically are 

Dtp = -p{V-v) (6) 
p {Dtv) = - 4Vp + vVv + + I) V ( V • t;) + (7) 

p (Dtc) =p [xV^c + V • (cVs)] + fc, (8) 

where all model parameters are constants. 

We note that none of the simplifying approximations we make above are necessary in principle. 
At the same time, not making such approximations requires knowing a number of physical prop- 
erties of the fluids, for example, the concentration dependence of the Soret coefficient St- Such 
information is difficult to obtain experimentally, and in any case, the known dependence is very 
weak and we believe it will not affect the results we present to within measurement or statistical 
error bars. Furthermore, accounting for the concentration dependence of the equation of state in 
the incompressible limit requires using variable-density low Mach number equations |48l HH] instead 
of the incompressible equations, since in general V • i; 7^ [50j. Extension of our algorithms to 



TpJS), although we emphasize that there are situations in 
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these variable-density variable-coefficient low Mach equations is possible but nontrivial, and will 
be considered in future work. 

II. FLUCTUATING HYDRODYNAMICS 

At mesoscopic scales the hydrodynamic behavior of fluids can be described with continuum 
stochastic PDEs of the Langevin type [HI [15j , as proposed by Landau and Lifshitz [16j and later 
justified by formal coarse-graining procedures [51]. Such equations can formally be justified as 
a central limit theorem for the Gaussian behavior of the thermal fluctuations around the mean, 
at least in certain simpler systems \52\ I53j. What emerges is that the (mesoscopic) thermal fluc- 
tuations can be described by the very same hydrodynamic equations describing the macroscopic 
behavior, linearized around the mean solution, and with added stochastic forcing terms that en- 
sure a fluctuation-dissipation balance principle [31]. Solving these linearized equations numerically 
requires first solving the deterministic equations for the mean, and then solving the fluctuating 
equations linearized around the mean. The linearization typically contains many more terms than 
the nonlinear deterministic terms due to the chain rule. Such a two step process is much more 
cumbersome then solving the nonlinear equations. Furthermore, a linearization has no hope of 
capturing any possible nonlinear feedback of the fluctuations on the mean flow, which is known to 
have physical significance pl|. 

Therefore, we follow an alternative approach in which the stochastic forcing terms are directly 
added to the nonlinear equations ( 6|7|8 ), but with an amplitude proportional to a parameter e that 



controls how far from linearity the equations are. In fluctuating hydrodynamics, to ensure mass 
and momentum conservation, the stochastic terms are the divergence of a stochastic flux, 

= e^V-S, and/e = eiv-*, (9) 

where the capital Greek letters denote stochastic fluxes that are modeled as white-noise Gaussian 
random fields. A detailed discussion of why there are no diffusive and stochastic fluxes in the density 
equation is given in Ref. [55j. For the linearized equations, we can fix e = 1, and the covariances 
of 5] and ^ can be derived from the fluctuation-dissipation balance principle, as explained well in 
the book [T7]. The covariance of the stochastic stress tensor Xl is not a positive definite matrix 
so there are many choices for how to express the stochastic stress, especially if additional bulk 
viscosity is included |56j . We have based our implementation on a formulation that requires the 
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fewest possible random numbers [511 [57] , 

S = + = + (^y^ - j Tr (W.) /, (10) 

* = ^/2xpM c(l-c) Wc (11) 

where W„ = (Wu + W^)/\/2 is a symmetric Gaussian random tensor field, and the \/2 in the 
denominator accounts for the reduction in variance due to the averaging. Here W,; and Wc 
are mutually-uncorrelated white-noise random Gaussian tensor and vector fields with uncorrelated 
components, 

{W\f{r, t)W^\r\ t')) = {6.,6,i) 6{t - mr - r') (12) 
(wf (r,t)wf (r',0) = {5^,)6{t-t')5{r-r'). (13) 

Similar covariance expressions apply in the Fourier domain as well if position r (time t) is replaced 
by wavevector k (wavefrequency uj), and {WaV^p) is replaced by (^WaW^y where star denotes 
complex conjugate (more generally, we denote an adjoint of a matrix or linear operator with a 
star). We recall that we take the temperature T to be spatially-constant. 

It is important to emphasize here that the non- linear fluctuating NS equations, with the white- 
noise stochastic forcing terms ([o]), are ill-defined because the solution should be a distribution 
rather than a function and the nonlinear terms cannot be interpreted in the sense of distributions. 
The nonlinear equations can be interpreted using a small-scale regularization (smoothing) of the 
stochastic forcing, along with a suitable renormalization of the transport coefficients [UiiSSj. Such 
a regularization is naturally provided by the discretization or coarse-graining [,57j length scale. As 
long as there are sufficiently many molecules per hydrodynamic cell the fluctuations will be small 
and the behavior of the nonlinear equations will closely follow that of the linearized equations 
of fluctuating hydrodynamics, which can be given a precise meaning |59j . This can be checked 
by reducing e to the point where the observed spatio-temporal correlations of the fluctuations 
begin to scale linearly with e, indicating nonlinear effects are negligible. In all of the simulations 
reported here, we have used e = 1 but have checked that using a very small e and then rescaling 
the covariance of the fluctuations by e^^ gives indistinguishable results to within statistical errors. 

Note that for the linearized equations the noise is additive since the covariance of the stochastic 
forcing terms is to be evaluated at the mean around which the linearization is performed. That 
is, in the linearized equations (e — )• 0) one should read ( [ll| as = \j2xpM c(l — c) Wc where c 
is the solution of the deterministic equations. Therefore, there is no Ito-Stratonovich difficulty in 
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interpreting the stochastic terms, and we use the (ambiguous) "Langevin" notation that is standard 
in the physics Uterature, instead of the differential notation more common in the hterature on 
stochastic differential equations. 

is a very challenging system of 

multiscale equations. Even a single stochastic advection-diffusion equation such as ([T]) is inher- 
ently multiscale because thermal fluctuations span the whole range of spatio-temporal scales from 
the microscopic to the macroscopic, specifically, all modes of the spatial discretization have a 
non-trivial stochastic dynamics that must be reproduced by the numerical method. Including the 
density equation ^ in the system of equations leads to fast sound wave modes that make the 
compressible equations stiff even in the deterministic setting. Finally, in most applications of inter- 
est the concentration diffusion is much slower than the momentum diffusion, leading to additional 
stiffness and multiscale nature of the equations, as we discuss further in Section |V) 



It is important to observe that even when linearized, (6|7p 



A. Incompressible Equations 



If density variations are negligible, p = po = const., we obtain the incompressible approximation 
to the hydrodynamic equations (6 |7|8 ) |17j . 



d^v = -V-K - V ■ {vv'^) + uV'^v + p'^fy 



(14) 
(15) 



where v = r///>, and v • Vc = V • (cv) and v • Vi> = V • [yv^^ because of incompressibility, 
V • t» = 0. Here 7' is the orthogonal projection onto the space of divergence-free velocity fields, 
"P = I — Q (T>Q ) ^ ^ X> in real space, where X>n = V • □ denotes the divergence operator and ^ = V 
the gradient operator. With periodic boundaries we can express all operators in Fourier space and 
= I — k^^{kk*) , where k is the wavenumber. 

Fluctuations can be included in the incompressible equations via the stochastic forcing terms 



p-^f^ =p-ieiV-*= V • \y/2exp-^M c(l - c) Wc 



Note that it is not necessary to include the stochastic pressure fluctuations Sp in (10) in the in 



compressible velocity equation (14) since the projection eliminates any non-zero trace component 



of the stress tensor. In our formulation, we use a strictly symmetric stochastic stress tensor in 
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the incompressible equations. This is based on physical arguments about local angular momentum 
conservation [60^ I61j. At the same time, the only thing that matters in the Fokker-Planck descrip- 
tion is the covariance of the stochastic forcing in the velocity equation T'D^s- This covariance is 
determined from the fluctuation dissipation balance principle, 

{{VVUs) {WT^sT) = VT> -DV = VCV, (16) 



where C is the vector Laplacian operator. Because and 2? have non-trivial null spaces, (16) 
does not uniquely determine the covariance of the stochastic stress. In fact, it can easily be shown 
by going to the Fourier domain that one can have a nonsymmetric component to the stochastic 



stress without violating (16). We believe that the stress tensor should be symmetric since we do 
not include an additional equation for the intrinsic spin (angular momentum) density. This is ap- 
propriate for fluids composed of "point" particles; however, recent molecular dynamics simulations 
have shown that for molecular liquids there can be non-trivial coupling between the linear and spin 
momentum densities |6lJ. While spin density has been included in the fluctuating hydrodynam- 
ics equations [60J at the theoretical level, we are not aware of any numerical simulations of such 
equations and do not consider an angular momentum equation in this work. 



B. Steady-State Covariances 

The means and spatio-temporal covariances of the fluctuating fields fully characterize the Gaus- 
sian solution of the linearized equations [28]. Of particular importance is the steady-state covariance 
of the fluctuating fields, which can be obtained for periodic systems by linearizing the equations in 
the fluctuations and using a spatial Fourier transform to decouple the different modes (wavevectors 
k). This steady-state covariance in Fourier space is usually referred to as a static structure factor 
in the physical literature, and represents the covariance matrix of the Fourier spectra of a typical 
snapshot of the fluctuating fields. Note that it is in principle possible to calculate the covariance 
of the fluctuations in non-periodic domains as well |62j : however, these tedious calculations offer 
little additional physical insight over the simple results presented below. We will present numeri- 
cal algorithms that can solve the fluctuating equations with non-periodic of boundary conditions; 
however, periodic conditions will be used to test the accuracy of the spatio-temporal discretization 
by comparing to the simple theory. 

At thermodynamic equilibrium, the fluctuations of the different hydrodynamic variables are 
uncorrelated and white in space, that is, the equilibrium variance is independent of the wavevector 
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k [28j, in agreement with equilibrium statistical mechanics \16\ I63j. Consider first the fluctuating 
isothermal compressible NS equations ( 6|7p ) linearized around a uniform steady-state, {p,v,c) = 
(Pq + 5p, vq + 5v, Co + Sc), T = Tq. Because of Galilean invariance, the advective terms vq ■ V (□) 
due to the presence of a background flow do not affect the equilibrium covariances (structure 
factors), which are found to be [17. ,28j 



Sr.r. 



6p) [^Sp 
6c] ( 5c 



PoksTo/cT 
PQ^kBToI 
Mp^^co{l-co). 



(17) 



At equilibrium, there are no cross-correlations between the different variables, for example, S^^v — 
(^{5c){5v)*'^ = 0. The equilibrium variance of the spatial average of a given variable over a cell of 
volume AV can be obtained by dividing the corresponding structure factor by AV, for example, the 
variance of the concentration is ^((5/9)^^ = poksTo/ (c|iAy). In the incompressible limit, ct — ?• oo, 
the density fluctuations vanish and p ^ pQ. 

Out of thermodynamic equilibrium, there may appear long-ranged correlations between the 
different hydrodynamic variables [17J. As a prototypical example of such non-equilibrium fluctua- 
tions, we focus on the incompressible equations ( 14|15 ) in the presence of an imposed concentration 
gradient Vc. The spatial non-uniformity of the mean concentration when there is a gradient breaks 
the translational symmetry and the Fourier transform no longer diagonalizes the equations. We 
focus our analysis and test our numerical schemes on a periodic approximation in which we linearize 
around a uniform background state {v, c) = {5v, cq + 5c), as suggested and justified in the physics 
literature on long-range nonequilibrium correlations \17\ \38\ [39l H5] . In such a periodic approxima- 
tion we cannot have a gradient in the steady-state average concentration cq but we can mimic the 
effect of the advective term v ■ Vcq with an additional term v ■ ( Vc) in the concentration equation. 
This is justified if the concentration gradient is weak, and leads to the linearized equations in a 
periodic domain. 



dti5v) = V 



{5v) + V • ( J2i^pQ^kBTo >V„ 



dt{5c) = -{Vc)-{5v)+xV^{5c) + V 



2xPo^Mco{l-co)Wc 



(18) 



In the Fourier domain (18) is a collection of stochastic differential equations, one system of linear 
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additive- noise equations per wavevector k, written in differential notation as 
d(6v] = -uk"^ (dv) dt + iJ2upQ^kBToVk- (dBi^^' 



d = - (Vc) • {dv^ dt - xk^ (Jcj dt + 2xPq^ M co{l - cq) k ■ [dSf'^^ , (19) 



where we used that 7' — Here Sv{t) is a tensor, and Bc{t) is a vector, whose components 
are independent Wiener processes. Note that the velocity equation is not affected by the concen- 



tration gradient. Given the model equations (19), the explicit solution for the matrix of static 
structure factors (covariance matrix) 

can be obtained as the solution of a linear system resulting from the stationarity condition dS = 0. 



For a derivation, see Eq. (30) in j^28j or Eq. (3.10) in [64j and also Eq. (41 ); below we simply quote 
the results of these straightforward calculations. 



1. Incompressible Velocity Fluctuations 

By considering the stationarity condition dSi,,v = it can easily be seen that the equilibrium 
covariance of the velocities is proportional to the projection operator, 



Sv,v = Po^kBToV = Po'kBTo [l - k-\kk*)] , 



(20) 



independent of the concentration gradient. In particular, the amplitude of the velocity fluctuations 
at each wavenumber is constant and reduced by one in comparison to the compressible equations, 

TV 5,,, = (^(6iy(fc)) = {d-l) Po^keTo, (21) 

where d is the spatial dimension. This is a reflection of the fact that one degree of freedom (i.e., 
one kBT/2) is subtracted from the kinetic energy due to the incompressibility constraint, which 



eliminates the sound mode. In Appendix |Bj we obtain the generalization of (20) to non-periodic 
systems, 



{vv*) = Po^kBTo {AV-^V) 



(22) 



An alternative way of expressing the result (22) is that all divergence- free modes have the 



same spectral power at equilibrium. That is, if the fluctuating velocities are expressed in any 
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orthonormal basis for the space of velocities that satisfy V • = 0, at equiUbrium the resulting 
random coefficients should be uncorrelated and have unit variance. This will be useful in Section 



IV A for examining the weak accuracy of the spatio-temporal discretizations of the incompressible 
equations. For periodic boundary conditions, such an orthonormal basis is simple to construct in 
the Fourier domain and a Fourier transform can be used project the velocity field onto this basis. 
In particular, for all wavevectors the projection of the velocity fluctuations onto the longitudinal 
mode 

^ = k \kx^ ky, kz] , (^^) 

where k = {k^ + ky + k'^Y^'^, should be identically zero, 

vi = (5v) • V = ^5vx + -r-Svy + ^6vz = k^^ (k • v) = 0. 
k k " k 

A basis for the incompressible periodic velocity fields can be constructed from the two vortical 
modes 

i>(2) = [kl+k^-^'\-ky,kx,Q], (24) 

^ = k (^k^ + kyj ^kxkz, kykz, — (^k^ + kyj'^ , (25) 

and the projection of the fluctuating velocities onto these modes has the equilibrium covariance 

(^2*2) = (^^3*3> = Po^^bTo, while {V2vl) = 0. (26) 

In two dimensions only v^^"^ and i)^'^^ are present, and i)^'^^ is the z component of the vorticity and 
spans the subspace of diverence-free velocities. The fact that the {d — 1) vortical modes have equal 



power leads to the velocity variance (21) 



2. Nonequilibrium Fluctuations 
When a macroscopic concentration gradient is present, the velocity fluctuations affect the con- 



centration via the linearized advective term (Vc) • v. Solving (19) shows an enhancement of the 



concentration fluctuations [651 proportional to the square of the applied gradient, 

Sc,c = Mp^' co(l - Co) + p^^^'f^^f^, i^^^' ^) (Vc)' , (27) 

where 9 is the angle between k and Vc, sin^ 9 = k'j_/k'^. Furthermore, there appear long-range 
correlations between the concentration fluctuations and the fluctuations of velocity parallel to the 
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concentration gradient, proportional to the applied gradient |1H I65j. 

5,,.,, = ((<rc)(^;)) = -J^^^^ (sin^ 0) V-c. (28) 

The power-law divergence for small k indicates long-range correlations between Sc and 6v and is 
the cause of the giant fluctuation phenomenon studied in Section [Vj 



III. SPATIO-TEMPORAL DISCRETIZATION 



Designing temporal discretizations for fluid dynamics is challenging even without including ther- 
mal fluctuations. When there is no stochastic forcing, our schemes revert to standard second-order 
discretizations and can be analyzed with existing numerical analysis techniques. Here we tackle 
the additional goal of constructing discretizations that, in a weak sense, accurately reproduce the 
statistics of the continuum fluctuations for the linearized equations. Note that achieving second- 
order weak accuracy is much simpler for linear additive-noise equations since in the linear case the 
solution is fully characterized by the means and the correlation functions (time-dependent covari- 
ances). In fact, one can use any method that is second-order in time in the deterministic setting 
and also reproduces the correct static (equal-time) covariance to second order, as explained in more 
detail in Ref. [28] . The deterministic order of accuracy can be analyzed using standard techniques, 
and the accuracy of the static covariances can be analyzed using the techniques described in Ref. 
|28j . We emphasize that the temporal integrators are only higher-order accurate in a weak sense 
for the linearized equations of fluctuating hydrodynamics \28\ H6] . 

Thermal fluctuations are added to a deterministic scheme as an additional forcing term that 
represents the temporal average of a stochastic forcing term over the time interval At and over the 
spatial cells of volume AV [28]. Because W is white in space and time, the averaging adds an addi- 
tional prefactor of {AV At) ' in front of the stochastic forcing. In the actual numerical schemes, 
a "realization" of a white-noise fleld W is represented by a collection W of normally-distributed 



random numbers with mean zero and covariance given by (13) or (12), with the identiflcation 

w ^ {AV Aty^i'^w. 



Speciflcally, the stochastic fluxes ( 10 ) are discretized as 



A realization of W is sampled using a pseudo-random number stream. The temporal dis- 
cretization of the stochastic forcing corresponds to the choice of how many realizations of W are 
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generated per time step, and how each reahzation is associated to specific points in time inside a 
time step (e.g., the beginning, mid-point, or end-point of a time step). The spatial discretization 
corresponds to the choice of how many normal variates to generate per spatial cell, and how to 
associate them with elements of the spatial discretization (e.g., cell centers, nodes, faces, edges). 
Once these choices are made, it is simple to add the stochastic forcing to an existing deterministic 
algorithm or code, while still accounting for the fact that white-noise is not like a classical smooth 
forcing and cannot be evaluated pointwise. 



As a first step in designing a spatio-temporal discretization for the compressible and incompress- 
ible equations of fluctuating hydrodynamics, we focus on the temporal discretization. We assume 
that the time step is fixed at At. The time step index is denoted with a superscript, for example, 
c" denotes concentration at time nAt and W"' denotes a realization of W generated at time step 
n. 

In the next section, we will describe our staggered spatial discretization of the crucial differential 
operators, denoted here rather generically with a letter symbol in order to distinguish them from 
the corresponding continuum operators. Specifically, let G be the gradient (scalar— )• vector), D 
the divergence (vector— T-scalar) , and L = DG the Laplacian (scalar— T-scalar) operator. When the 
divergence operator acts on a tensor field F such as a stress tensor cr, it is understood to act 
component-wise on the x, y and z components of the tensor. Similarly, the gradient and Laplacian 
act component-wise on a vector. An important property of the discrete operators that we require 
to hold is that the divergence operator is the negative adjoint of the gradient, D = —G*. This 
ensures that the scheme satisfies a discrete version of the continuous property, 



for any scalar field w{r). 

We define the weak order of accuracy of a temporal discretization in terms of the mismatch 
between the steady-state covariance of the continuum and the discrete formulations. With peri- 
odic boundary conditions this would be the mismatch between the Fourier spectrum of a typical 
snapshot of the true solution and the steady-state discrete spectrum of the numerical solution |28j . 
This mismatch is typically of the form 0{l^t^) for some integer A; > 1, implying that for sufficiently 
small time steps the discrete formulation reproduces the steady-state covariance of the continuum 



A. Temporal Discretization 





V ■ Vu; dr if 17 • uqq = or t; is periodic 
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formulation. Note that for the Unearized equations a certain order of deterministic temporal accu- 
racy, combined with equal or higher order of accuracy of the steady-state covariances, implies the 
same order of accuracy for all temporal correlations. A theoretical analysis of the weak accuracy of 
the temporal discretizations used in this work can be performed using the tools described in Ref. 
|28j with some straightforward extensions [1^; here we simply state the main results and verify 
the order of weak accuracy numerically. 



1. Compressible Equations 

Denoting the fluctuating field with Q = {p,v,c), the fluctuating compressible NS equations 
( 6|7|8 ) can be written as a general stochastic conservation law, 



dtQ = -D [F{Q- 1) - Z{Q, W)] , (30) 
where D is the divergence operator (acting component- wise on each flux), F{Q;t) is the deter- 



ministic and = [0, S, ^'J is the discretization of the stochastic flux (29). We recall that the 
stochastic forcing amplitude is written as multiplicative in the state; however, in the linearized 
limit of weak fluctuations the strength of the stochastic forcing only depends on the mean state, 
which is well-approximated by the instantaneous state, Q{t) ^ Q{t). Following we base our 



temporal discretization of (30) on the (optimal) three-stage low-storage strong stability preserving 
|66j (originally called total variation diminishing |67| ) Runge-Kutta (RK3) scheme of Gottlieb and 
Shu, ensuring stability in the inviscid limit without requiring slope-limiting. The stochastic terms 
are discretized using two random fluxes per time step, as proposed in Ref. [28j. This discretization 
achieves third-order weak accuracy [46j for linear additive-noise equations, while only requiring the 
generation of two Gaussian random flelds per time step. 

For each stage of our third-order Runge-Kutta scheme, a conservative increment is calculated 

as 

AQ(Q, W; t) = -At DF{Q; t) + At DZ{Q, W). 

Each time step of the RK3 algorithm is composed of three stages, the first one estimating Q at time 
t = {n + l)At, the second at t = (n + ^)At, and the final stage obtaining a third-order accurate 
estimate at t = (n + l)At. Each stage consists of an Euler-Maruyama step followed by a weighted 
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averaging with the value from the previous stage, 

;n+l 



Q 



=g" + AQ(Q",Vr^ nAt) 
^-Q" + I [Q +AQ(^Q^,W^;{n + l)At 



(31) 



where the stochastic fluxes between diff'erent stages are related to each other via 



(32) 



and 'W\ and VF^ are two independent realizations of W that are generated independently at each 
RK3 step. In this work we used the weights derived in Ref. |28j based on a linearized analysis, 
w\ = — "v/S, W2 = \/3 and = 0. More recent analysis based on the work in |68j shows that second- 
order weak accuracy is achieved for additive-noise nonlinear stochastic differential equations using 
the weights jl6] 



Wl 



{2V2-V3) 



W2 



-4^/2-3^/3) (^/2 2 V3) 

— , and ^3 = - 



5 5 10 

For the types of problems studied here nonlinearities play a minimal role and either choice of the 
weights is appropriate. 



2. Incompressible Equations 



The spatially-discretized equations ( 14|15 ) can be written in the form 

dtV + Gir = Av{v,c) + uLv + p^^ f^, 
dtc = Ac{v,c) + xLc + p~'^fc, 
s.t. Dv = 0, 



where A{v,c) represent the non-diffusive deterministic terms, such as the advective and Soret 
forcing terms , as well as any additional terms arising from gravity or other effects. Fluctuations 
are accounted for via the stochastic forcing terms 



D 



2evp ^ksT 
AV At 



and p Vc (c, Wc 



D 



2exp-^M c(l - c) 
Ay At 
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We base our temporal discretization on the second-order semi-implicit deterministic scheme of 
Griffith [3l] , a predictor-corrector method in which the predictor step combines the Crank-Nicolson 
method for the diffusive terms with the Euler method for the remaining terms, 



At 



At 

s.t. Dv^^^^ = 0. (33) 

The corrector stage combines Crank-Nicolson for the diffusive terms with an explicit second-order 
approximation for the remaining deterministic terms, 

I'^+l _ -J," , „,i /^,n+l, n\ 1 



+ G7r"+2 = A^'^ +uL\ + p~'f. 



= Ac ^ +xL[ ]+ P fc ' 



s.t. Dv''+'^ = 0. (34) 

Unlike a fractional-step scheme that splits the velocity and pressure updates [Ml EO] , this approach 
simultaneously solves for the velocity and pressure and avoids the need to determine appropriate 
"intermediate" boundary conditions. Importantly, no spurious boundary modes [711 [72] arise due 
to the implicit velocity treatment even in the presence of physical boundaries, which is especially 
important for fluctuating hydrodynamics since all of the modes are stochastically forced |46| . 



The concentration equation in (34) [and similarly in (33)] is a linear system for c""*"^ that ap- 
pears in standard semi-implicit discretizations of diffusion and is solved using a standard multigrid 



method. The velocity equation in ( |34[ ) [and similarly in (33)] is a much harder "saddle-point" 
systems of linear equations to be solved for the variables 1?"+^ and vr"^?. This time-dependent 
Stokes problem is solved using a Krylov iterative solver as described in detail in Ref. [34j- The 
ill-conditioning of the Stokes system is mitigated by using a projection method (an inhomogeneous 
Helmholtz solve for velocity followed by a Poisson solve for the pressure) as a preconditioner. With 
periodic boundary conditions solving the Stokes system is equivalent to a projection method, that 
is, to an unconstrained step for the velocities followed by an application of the projection operator. 
If physical boundaries are present, then the projection method is only an approximate solver for 
the incompressible Stokes equations; however, the "splitting error" incurred by the approximations 
inherent in the projection method is corrected by the Krylov solver. 

The nonlinear terms are approximated in the corrector stage using an explicit trapezoidal rule, 

aT'^ = ^ [A„(^", c") + A„(i3"+i, 5"+!)] , (35) 
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which is the (optimal) two-stage strong stabihty preserving Runge-Kutta method [66], and is thus 
generally preferable for hyperbolic conservation laws. For the stochastic forcing terms, we employ 
a temporal discretization that uses one random flux per time step, 



fv' = /. [W,) and /r ^ = (^c"+2, 1^"^ 



where cp'~^2 = (c" + c"^"^) /2 but we again emphasize that the dependence of fc ^ on the instan- 
taneous state is not important in the weak- noise (linearized) setting. It can be shown that 
this temporal discretization is second-order weakly accurate for additive-noise nonlinear stochastic 
differential equations [46j. More importantly, the Crank-Nicolson method balances the numerical 
dissipation with the stochastic forcing identically in the linear setting. This important property 
allows our time stepping to under-resolve the fast dynamics of the small-wavelength fluctuations 
while still maintaining the correct spectrum for the fluctuations at all scales. While this surprising 
fact has already been verified (in a simplified setting) in the Appendix of Ref. [28] and also in Ref. 
|73j . we give a different derivation in Appendix [A} A more detailed analysis will be presented in a 
forthcoming paper 



The linearized equations ( 18 ) have additional structure that enables us to simplify the predictor- 
algorithm. Firstly, the momentum equation is independent of the concentration equation(s), 
A„(u,c) = 0, and the corrector step of the velocity equation is redundant since it simply re- 
peats the predictor step, i)"^^ = Therefore, we only need to do one Stokes solve per time 
step. Furthermore, only velocity enters in the linearized concentration equation, Ac{v,c) = Av, 
and therefore 

Ac ^ = Av"-^2 = A- = Av"-^2 

2 

can be calculated without performing a predictor step for the concentration. This variation of the 
time stepping is twice as efficient and can be thought of a split algorithm in which we first do a 
Crank-Nicolson step for the velocity equation, 

— + G7r"+2 = „L ( ^ j + p-^f^ [W^j such that £>i;"+i = 0, (36) 

and then a Crank-Nicolson step for the concentration equation using the midpoint velocity to 
calculate advective fluxes. 



-J + xL J + p-'fc {c\ W:) . (37) 



At V 2 

Because of the special structure of the equations, the split algorithm is equivalent to the tradi- 
tional Crank-Nicolson method applied to the coupled velocity-concentration system in which both 
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advection and diffusion are treated semi-implicitly. This observation, together with the derivation 
in Appendix [Aj shows that our temporal discretization gives the correct steady-state covariances 
for any time step size At, although it does not reproduce the correct dynamics for large At. This 
property will prove very useful for the simulations of giant fluctuations reported in Section |V| 



B. Spatial Discretization 

We now consider spatial discretization of the equations of fluctuating hydrodynamics on a regu- 
lar Cartesian grid, focusing on two dimensions for notational simplicity. The spatial discretization 
is to be interpreted in the finite-volume sense, that is, the value of a fluctuating field at the center 
of a spatial cell of volume AV represents the average value of the fluctuating field over the cell. We 
explicitly enforce strict local conservation by using a conservative discretization of the divergence. 
Specifically, the change of the average value inside a cell can always be expressed as a sum of fluxes 
through each of the faces of the cell, even if we do not explicitly write it in that form. 

Consider at first a simplified form of the stochastic advection-diffusion equation for a scalar 
concentration field 



dtc 



(38) 



where v{r,t) is a given advection velocity. We note that for incompressible fiow, we can split the 
stochastic stress tensor into a vector Wx corresponding to the flux for v^, and a vector Wy 
corresponding to Vy. We can then view the velocity equation as a constrained pair of stochastic 
advection-diffusion equations of the form (38), one equation for Vx and another for Vy. We will 



discuss the generalization to compressible flow in Section III B 5 

The spatial discretization described in this section is to be combined with a suitable stable 
temporal discretization, speciflcally, the temporal discretization that we employ was described in 



Section [III A We consider here the limit of small time steps. At — )• 0, corresponding formally to a 



semi-discrete "method of lines" spatial discretization of the form 

dc 



D 

dt 



{-Uc + xGc) + V2x/ (AF At)Wc , (39) 



where c = {qj} is a flnite- volume representation of the random field c{r,t). Here, D is a. con- 
servative discrete divergence, G is a discrete gradient, and U = U (v) denotes a discretization of 
advection by the spatially-discrete velocity field v, and Wc denotes a vector of normal variates 
with specifed covariance Cw = i'^cW*). 
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1. Discrete Fluctuation- Dissipation Balance 

We judge the weak accuracy of the spatial discretization by comparing the steady-state co- 
variance of the spatially-discrete fields to the theoretical covariance of the continuum fields in the 
limit — 7- [28j. Ignoring for a moment constraints such as incompressibility, at thermodynamic 
equilibrium the variance of the discrete fields should be inversely proportional to Ay and values 
in distinct cells should be uncorrelated 

Ce = (cc*> = Sc,c (Ay-1/) . (40) 

For periodic systems this means that the spectral power of each discrete Fourier mode be equal to 



the continuum structure factor, 5c,c = 1 for the model equation (38) [see also (17)], independent 
of the wavenumber. 

A spatial discretization that gives the correct equilibrium discrete covariance is said to satisfy 
the discrete fluctuation-dissipation balance (DFDB) condition |28| I46j. The condition guarantees 
that for sufficiently small time steps the statistics of the discrete fluctuations are consistent with the 
continuum formulation. For larger time steps, the difference between the discrete and continuum 
covariance will depend on the order of weak accuracy of the temporal discretization |74j . A simple 
way to obtain the DFDB condition is from the time stationarity of the covariance. For the model 



equation (38) we obtain the linear system of equations for the matrix Cc 



dC d (rr*\ 

= = D{-U + xG) Ce + Ce [D {-U + xG)r + 2xAV~'DCwD^ = 0, (41) 



whose solution we would like to be given by (40), specifically, Cc = Ay I. Considering first the 



case of no advection, U = 0, we obtain the DFDB condition 

DG + {DCy = -2DCwD*- (42) 
Consider first the case of periodic boundary conditions. A straightforward way to ensure the 



condition ( 42 ) is to take the components of the random flux Wc to be uncorrelated normal variates 
with mean zero and unit variance, Cw = aud also choose the discrete divergence and gradient 
operators to be negative adjoints of each other, G = —D*, just as the continuum operators 



are [251 1211 EH (see Eq. 43). Alternative approaches and the advantages of the above "random 



flux" approach are discussed in Ref. jG^. As we will demonstrate numerically in Section IV, the 
staggered discretization of the dissipative and stochastic terms described below satisfles the discrete 
fluctuation-dissipation balance for both compressible and incompressible flow. 
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In the continuum equation (38), the advective term does not affect the fluctuation-dissipation 
balance at equilibrium; advection simply transports fluctuations without dissipating or amplifying 
them. This follows from the skew-adjoint property 

/ If [V • {cv)] dr = — c [V • {wv)] dr if V • = and v ■ ngn = or t) is periodic, (43) 
Jn Jn 

which holds for any scalar field w{r). In particular, choosing w = c shows that for an advection 

equation dtc = —V • (cv) the "energy" / c^dr/2 is a conserved quantity. To ensure that the 



discrete fluctuation-dissipation balance (41) is satisfied, the matrix DUCc, or more precisely, the 
discrete advection operator S = DU should be skew-adjoint, = —S. Specifically, denoting with 
c - w = J2ij the discrete dot product, we require that for all w 

w ■ [{DU) c] = -c • [{DU) w] (44) 

if the advection velocities are discretely-divergence free, {DU) 1 = 0, where 1 denotes a vector 
of all ones. Note that this last condition, 51 = 0, ensures the desirable property that the advec- 
tion is constant-preserving, that is, advection by the random velocities does not affect a constant 
concentration field. 

For incompressible fiow, the additional constraint on the velocity Dv = needs to be taken 



into account when considering discrete fiuctuation-dissipation balance. In agreement with (22), we 
require that the equilibrium covariance of the discrete velocities be 

{vv*) = Po'kBTo (Ay-ip) , (45) 

where P is the discrete projection operator 

F = I-G {DG)'^ D = I- D* {DD*)''^ D. 



With periodic boundary conditions, (45) implies that the discrete structure factor for velocity is 
"S*,,,,; = Pq^UbTq^. In particular, the variance of the velocity in each cell is in agreement with 
the continuum result, since TrP = Tr7' = (i— 1. More generally, for non-periodic or non-uniform 
systems, we require that for sufficiently small time steps all discretely-incompressible velocity modes 



are equally strong at equilibrium [16]. In Appendix|B|we generalize the DFDB condition (42 ) to the 
incompressible (constrained) velocity equation, and show that there are no additional conditions 
required from the discrete operators other than the duality condition on the divergence and gradient 
operators, G = —D*. 
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2. Staggered Grid 



A cell-centered discretization that is of the form (39) and satisfies the discrete fluctuation- 
dissipation balance (DFDB) condition was developed for compressible flow in Ref. [2H]. Extending 
this scheme to incompressible flow is, however, nontrivial. In particular, imposing a strict discrete 
divergence-free condition on collocated velocities has proven to be difficult and is often enforced only 



approximately [73], which is inconsistent with (45), as we explain in Appendix pj An alternative 



is to use a staggered grid or "MAC" discretization, as first employed in projection algorithms for 
incompressible flow [76j. In this discretization, scalars are discretized at cell centers, i.e., placed 
at points (i, j), while vectors (notably velocities) are discretized on faces of the grid, placing the x 
component at points (i+1/2, j), and the y component at (i, j + 1/2). Such a staggered discretization 
is used for the fluxes in Ref. [28|, the main difference here being that velocities are also staggered. 

In the staggered discretization, the divergence operator maps from vectors to scalars in a locally- 
conservative manner, 

V • ^ ^ {DvV , = Ax-i ( . - + Ay-^ ( v^"^ , - v^''^ , \ . 

The discrete gradient maps from scalars to vectors, for example, for the x component: 



(^) _ A^-l 

2 ' 



It is not hard to show that with periodic boundary conditions G = —D* as desired. The resulting 
Laplacian L = DG is the usual 5-point Laplacian, 

which is negative definite except for the expected trivial translational zero modes. The velocities 
Vx and Vy can be handled analogously. For example, Vx is represented on its own finite- volume 
grid, shifted from the concentration (scalar) grid by one half cell along the x axis. The divergence 
D^^\ gradient G^^^ and Laplacian L^^^ are the same MAC operators as for concentration, but 
shifted to the x-velocity grid. 

For the compressible equations, there is an additional dissipative term in (|4]) that involves 
V(V-i7). This term is discretized as written, GDv, which can alternatively be expressed in 
conservative form. When viscosity is spatially-dependent, the term V- (r/Vi^) should be discretized 
by calculating a viscous flux on each face of the staggered grids, interpolating viscosity as needed 
and using the obvious second-order centered differences for each of the terms dxVx, dxVy, dyVy and 
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dyVx- For a collocated velocity grid the mixed derivatives dxVy and dyVx, and the corresponding 
stochastic forcing terms, do not have an obvious face-centered discretization and require a separate 
treatment [2H] . 



3. Stochastic Fluxes 



The stochastic flux Wc, like other vectors, is represented on the faces of the grid, that is, Wc is 
a vector of i.i.d. numbers, one number for each face of the grid. To calculate the state-dependent 



factor c(l — c) that appears in (29) on the faces of the grid, concentration is interpolated from 
the cell centers to the faces of the grid. At present, lacking any theoretical analysis, we use a simple 



arithmetic average (47) for this purpose. 

The stochastic momentum flux Wv is represented on the faces of the shifted velocity grids, which 
for a uniform grid corresponds to the cell centers and the nodes (z+ |, j + ^) of the grid PU] . 

(x) (y) 

Two random numbers need to be generated for each cell center, W"- • and W"- ^ , corresponding to 

the diagonal of the stochastic stress tensor. Two additional random numbers need to be generated 

(x) (v) 

for each node of the grid, W: { . i and W \ . i , corresponding to the off-diagonal components. 

^+2^+2 *+2'-J"^2 

In three dimensions, the three diagonal components of the stochastic stress are represented at the 
cell centers, while the six off-diagonal components are represented at the edges of the grid, two 

(x) (v) 

random numbers per edge, for example, W. \ i and W.^\ • , i i.- 

For the incompressible equations one can simply generate the different components of as 
uncorrelated normal variates with mean zero and unit variance, and obtain the correct equilibrium 
covariances. Alternatively, each realization of the stochastic stress can be made strictly symmetric 



and traceless as for compressible flow, as specified in (10). Because of the symmetry, in practice 
for each node or edge of the grid we generate only a single unit normal variate representing the two 
diagonally-symmetric components. For each cell center, we represent the diagonal components by 
generating d independent normal random numbers of variance 2 and then subtracting their average 
from each number. Note that for collocated velocities a different approach is required because the 
diagonal and diagonally-symmetric components of the stress tensor are not discretized on the same 
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4- Advection 

We now consider skew-adjoint discretizations of the advection operator S = DU on a staggered 
grid. This problem has been considered in a more general context for the purpose of constructing 
stable methods for turbulent flow in Ref. \n\ I78j : here we focus on a simple second-order cen- 
tered discretization. The importance of the skew-adjoint condition in turbulent flow simulation 
is that it leads to strict discrete energy conservation for inviscid flow, which not only endows the 
schemes with long-time stability properties, but also removes undesirable numerical dissipation. 
Conservation of the discrete kinetic energy Ek = p {v ■ v) /2 is also one of the crucial ingredients 
for fluctuation-dissipation balance, i.e., the requirement that the Gibbs-Boltzmann distribution 
exp [—El:/ {^^bT)] be the invariant distribution of the stochastic velocity dynamics [19 \ [25 l [79]. 

Consider first the spatial discretization of the advective term DUc in the concentration equa- 
tion. Since divergence acts on vectors, which are represented on the faces of the grid, Uc should 
be represented on the faces as well, that is, U is a linear operator that maps from cell centers to 
faces, and is a consistent discretization of the advective flux cv. If we deflne an advection velocity 
u on the faces of the grid, and also define a concentration c on each face of the grid, then the 
advective fiux can directly be calculated on each face. For example, for the x faces: 

^ iUc)^y = (46) 

For concentration we can take u = v, since the velocity is already represented on the faces of the 
scalar grid. Simple averaging can be used to interpolate scalars from cells to faces, for example, 

Cj+ij = ^(ci+ij + Cij), (47) 

although higher-order centered interpolations can also be used p8|. 



As discussed in Section IIIB 1 , we require that the advection operator be skew adjoint if DUl 



Du = 0. Our temporal discretization of the incompressible equations ( 33|34 ) ensures that a 



discretely divergence-free velocity is used for advecting all variables. The case of compressible flow 



will be discussed further in Section III B 5 In the incompressible case, S = DU can be viewed as 



a second-order discretization of the "skew-symmetric" form of advection [77] 

c 1 
v-Vc + -V -v = - [V ■ (cv) + V ■ Vc] . 



Namely, using ( 46 ) we obtain 
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and rewrite the x term using (47) as 



{x) {x) 
, 1 .C,- , 1 - U. \ .C,_l 

«+9J *+2J 2-9^ * 2'-? 



+ c 



(x) {x) 
U. I . — U. 1 . 

«+2J «-2J 



and similarly for the y term, to obtain 



[DUc] 



(Sc 



where is a centered discretization of [V • (cv) + v ■ Vc] /2, 



5c 



,(y) 



(48) 



• (49) 



w, and, assuming 



Since the advection velocity is discretely divergence free, S = S. 

It is not hard to show that S is skew-adjoint. Consider the x term in Sc 
periodic boundary conditions, shift the indexing from i to i — 1 in the first sum and from i to i + 1 
in the second sum, to obtain 



Or 

1 

1-2,3 '■' 



Therefore, S is skew-adjoint, ( Sc] • it) = — c • ( Sw ) . A similar transformation can be performed 



with slip or stick boundary conditions as well. These calculations show that (44) holds and thus 



the discrete advection operator is skew-adjoint, as desired. Note that the additional terms in 



(15) due to the Soret effect can be included by advecting concentration with the effective velocity 



■Wadv = v - xSt'VT. 

The same approach we outlined above for concentration can be used to advect the velocities as 
well. Each velocity component lives on its own staggered grid and advection velocities are needed 
on the faces of the shifted grid, which in two dimensions corresponds to the cell centers and the 
nodes of the grid. The velocity is advected using an advection velocity field u^^^ that is obtained 
via a second-order interpolation of v, 



u. 



(x) 



+ , 1 



1 

2 

2 " 



and similarly for the other components. It is not hard to verify that the advection velocity ii^^) is 
discretely divergence-free if v is: 

1 



showing that D^^^ii^^) = if Dv = 0. Therefore, the shifted advection operator S^^^ = D^^^U^^^ 
is also skew-adjoint, as desired. 
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5. Compressible Equations 

It is instructive at this point to summarize our spatial discretization of the incompressible 
equations ( 14|15 ), before turning to the compressible equations. The concentration equation (15) 



is discretized as 

dc 



-DUc + xDGc + D^, (50) 



where U is given by (46) with advection velocity u = v — xSt^T. For the x component of the 
velocity we use the spatial discretization 

and similarly for the other components, and the pressure ensures that Dv = 0. 

Our staggered spatial discretization of the compressible equations ( 6|7|8 ) is closely based on 



the discretization described above for the incompressible equations. An important difference is 
that for compressible flow we use the conservative form of the equations, that is, we use the mass 
density p, the momentum density j = pv and the partial mass density pi = cp as variables. The 
momentum densities are staggered with respect to the mass densities. Staggered velocities are 
defined by interpolating density from the cell centers to the faces of the grid, for example, 

which implies that Dj = DUp. 

The density equation ([6]) is discretized spatially as 

t = -OUp, (51) 

while for the concentration equation ([s]) we use 

^ = -DUpi + poXoDGc + (52) 

where we assume that px = PoXo is constant. For the x component of the momentum density we 
use 

^ = -D^-)U^^)j^ - 4 [Gp)^ + ryD(")G(")i),. + (C + ^) {GDv)^ + D^^^'E^^\ (53) 

and similarly for the other components. The spatio-temporal discretization ensures strict local 
conservation of p, j and pi. 
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The discretization ( 51|52|53 ) satisfies discrete- fluctuation dissipation balance at equilibrium, 
specifically, the equilibrium covariances of velocity and density are (vv*) = PQ^ksTo I and (pp*) = 



poksTo/cj^ I, in agreement with the continuum spectra given in (17). Linearizing the semi-discrete 



density equation (51 ) around an equilibrium state (p, v) = (pg + 6p, vq + 6v) with Dvq = gives 

d{5p) 



dt 



+ So (dp) = -po [D {6v)] . 



Recall that the operator Sq, defined by (49) with u = vq, is skew-adjoint, and the fluctuations in 
density are thus controlled by the coupling with the velocity fluctuations. For simplicity, consider 
this coupling for the case of a fluid at rest, vq = and thus Sj = po {Sv). Linearizing the momentum 



update (53) and focusing on the coupling with the density fluctuations, we obtain 

— — h advection = —Pq^c^ [G {6p)] + dissipation and forcing. 

Fluctuation-dissipation balance requires the skew-symmetry property Lp .^ (vv*) = — {pp*) L'^ p, 
where Lp_„ = —poD the operator in front of 5v in the density equation, and L^^p = —c^G is 
the operator in front of 5p in the velocity equation. This skew-symmetry requirement is satisfied 
because of the key duality property D = —G*. This demonstrates the importance of the duality 
between the discrete divergence and gradient operators, not just for a single advection-diffusion 
equation, but also for coupling between the different fluid variables. In future work, we will explore 
generalizations of the concept of skew- adjoint discrete advection to the nonlinear compressible 
equations [56l [78] . 



6. Boundary Conditions 

Non-periodic boundary conditions, specifically, Neumann or Dirichlet physical boundaries, can 
be incorporated into the spatial discretization by modifying the discrete divergence, gradient and 
Laplacian operators near a boundary. This needs to be done in a way that not only produces an 
accurate and robust deterministic scheme, but also ensures fluctuation-dissipation balance even in 
the presence of boundaries. Here we extend the approach first suggested in an Appendix in Ref. 
|13j to the staggered grid. It can be shown that the inclusion of the (discrete) incompressibility 
constraint does not affect the fluctuation-dissipation balance when an unsplit Stokes solver is 
employed in the temporal integrator }46j . 

We assume that the physical boundary is comprised of faces of the grid. Since only the direction 
perpendicular to the wall is affected, we focus on a one-dimensional system in which there is a 
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physical boundary between cells 1 and 0. For the component of velocity perpendicular to the wall, 
some of the grid points are on the physical boundary itself and those values are held fixed and 
not included as independent degrees of freedom. For the second-order spatial discretization that 
we employ no values in cells outside of the physical domain are required. Therefore, no special 
handling at the boundary is needed. 

For cell-centered quantities, such as concentration and components of the velocity parallel to 
the wall, the boundary is half a cell away from the cell center, that is, the boundary is staggered. 
In this case we use the same discrete operators near the boundaries as in the interior of the domain, 
using ghost cells extending beyond the boundaries to implement the finite-difference stencils near 
the boundaries. One can think of this as a modification of the stencil of the Laplacian operator 
L near boundaries, specifically, when boundaries are present the dissipative operator L ^ DG 



but rather L = DG, where G is a modified gradient. Repeating the calculation in (41) for the 
spatially-discretized model equation 

^ = xLc+ V2x/ iAVAt)DW. 



leads to a generalization of the DFDB condition (42), assuming L* = L, 

xLCa + xCcL* = 2xAV-^L = -2xAV-^ DCwD* ^ L = -DCwD*. (54) 

Consider first a Neumann condition on concentration, dc{0)/dx = 0. This means that a no-flux 
condition is imposed on the boundary, and therefore for consistency with physical conservation the 
stochastic flux on the boundary should also be set to zero, Wi = 0. The ghost cell value is set 

2 

equal to the value in the neighboring interior cell (reflection), cq = ci, leading to 

(DW)^ = Ax-^ Wa , (Gc) ^ = 0, (Lc)i = Ax~'^ (c2 - ci) . (55) 



2 



If we exclude points on the boundary from the domain of the divergence operator, which is also 
the range (image) of the gradient operator, then it is not hard to see that the duality condition 
D* = —G continues to hold. We can therefore continue to use uncorrelated unit normal variates 



for the stochastic fluxes not on the boundary, = / in (54). 

If a Dirichlet condition c(0) = is imposed, then the ghost cell value is obtained by a linear 
extrapolation of the value in the neighboring interior cell (inverse reflection), cq = — ci, leading to 

(DW)^ = Ax-^ (W3 - Wi^ , (GcJ ^ = Ax-^ (2ci) , [Lc]-^ = Ax'^ (cs - 3ci) . (56) 

The duality condition is no longer satisfied, D* ^ —G, but it is not hard to show that the 



fiuctuation-dissipation balance condition (54) can be satisfied by simply doubling the variance of 
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the stochastic flux on the boundary, {WiWt ) = 2. Note that the Laplacian (56) is not formany 



second-order accurate at the boundary, however, its normal modes (eigenvectors) can be shown to 
correspond exactly to the normal modes of the continuum Laplacian and have decay rates (eigen- 
modes) that are second-order accurate in Ax^, and in practice pointwise second-order accuracy is 
observed even next to the boundary. Formal second-order local accuracy can be obtained by using 
a quadratic extrapolation for the ghost cell, cq = —2ui + M2/3 and (Lc)-^ = Ax~^ (4c2/3 — 4ci), 
however, this requires a more complicated handling of the stochastic fluxes near the boundary as 
well. 

In summary, the only change required to accommodate physical boundaries is to set the variance 
of stochastic fluxes on a physical boundary to zero (at Neumann boundaries), or to twice that used 
for the interior faces (at Dirichlet boundaries). For density in compressible flows, the ghost cell 
values are generated so that the pressure in the ghost cells is equal to the pressure in the neigh- 
boring interior cell, which ensures that there is no unphysical pressure gradient in the momentum 
equation across the interface. There is also no stochastic mass flux through faces on the boundary 
independent of the type of boundary condition at the wall. For incompressible flow the gradient of 
pressure is discretized as Gn = —D*tt even in the presence of stick or slip boundary conditions for 
velocity; more complicated velocity-stress or open |27j boundary conditions are simple to handle 
(in principle) with the projection-preconditioner solvers. 

IV. IMPLEMENTATION AND NUMERICAL TESTS 

We now describe in more detail our implementations of the spatio-temporal discretizations 



described in Section III, and provide numerical evidence of their ability to reproduce the correct 
fluctuation spectrum in uniform flows with periodic boundary conditions. A less trivial application 
with non-periodic boundaries is studied in Section [V] 

We consider here a uniform periodic system in which there is a steady background (mean) 
flow of velocity vq. Unlike the continuum formulation, the discrete formulation is not Galilean- 
invariant under such uniform motion and the covariance of the discrete fluctuations is affected by 
the magnitude of vq. The stability and accuracy of the spatio-temporal discretization is controlled 
by the dimensionless CFL numbers 

VAt ^ i^At x^t 
where V = ct (isothermal speed of sound) for low Mach number compressible flow, and V = || I'd 1 1 00 
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for incompressible flow, and typically x ^- The explicit handling of the advective terms places a 
stability condition a < 1 , and the explicit handling of diffusion in the compressible flow case requires 
/3) /3c ^ 1/2"^, where d is the dimensionality. The strength of advection relative to dissipation is 
measured by the cell Reynolds number r = a/ (3 = V / [u/S.x). 

To characterize the weak accuracy of our methods we examine the discrete Fourier spectra 
of the fluctuating fields at equilibrium, and compare them to the continuum theory discussed in 



Section II B for all discrete wavenumbers k. We use subscripts to denote which pair of variables is 
considered, and normalize each covariance so that for self-correlations we report the relative error 
in the variance, and for cross-correlations we report the correlation coefficient between the two 
variables. For example, the non-dimensionalized static structure factor for concentration is 

where c{k) is the discrete Fourier transform of the concentration. Note that an additional factor 
equal to the total number of cells may be needed in the numerator depending on the exact definition 
used for the discrete Fourier transform [25]. Similarly, the cross-correlations between different 
variables need to be examined as well, such as for example, 

Sc,v = I {cv*) . 

^[Mpoico(l-co)] (po'^B^o) 

For staggered variables the shift between the corresponding grids should be taken into account as 
a phase shift in Fourier space, for example, exp(A;^Ax/2) for Vx- For a perfect scheme, S^c = 1 
and Sc,v = for all wavenumbers, and discrete fluctuation-dissipation balance in our discretization 
ensures this in the limit At — )■ 0. Our goal will be to quantify the deviations from "perfect" for 
several methods, as a function of the dimensionless numbers a and (3. 

A. Incompressible Solver 



We have implemented the incompressible scheme described in Sections III A 2 and III B using the 
IB AMR software framework [SOj, an open-source library for developing fluid-structure interaction 
models that use the immersed boundary method. The IBAMR framework uses SAMRAI [81] to 
manage Cartesian grids in parallel, and it uses PETSc [82j to provide iterative Krylov solvers. The 
majority of the computational effort in the incompressible solver is spent in the linear solver for 
the Stokes system; in particular, in the projection-based preconditioner, the application of which 
requires solving a linear Poisson system for the pressure, and a modifled linear Helmoltz system for 
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the velocities and the concentrations [33] . For small viscous CFL numbers /? ^ 1 the Poisson solver 
dominates the cost, however, for (3^1 the two linear systems become similarly ill-conditioned 
and require a good preconditioner themselves. We employ the hypre library [HS] to solve the linear 
systems efficiently using geometric multigrid solvers. 

For incompressible flow, one could directly compare the spectrum of the velocities {vv*) to the 



spectrum of the discrete projection operator P (see Section III B 1 ) . It is, however, simpler and more 
general to instead examine the equilibrium covariance of the discrete modes forming an orthonormal 
basis for the space of discretely divergence free modes. The amplitude of each mode should be unity 
for all wavenumbers, even if there are physical boundaries present, making it easy to judge the 
accuracy at different wavenumbers. For periodic boundary conditions a discretely-orthogonal basis 
is obtained by replacing the wavenumber k = [kx, ky, k^) in ( 23|24|25 ) by the effective wavenumber 
k that takes into account the centered discretization of the projection operator, for example, 

~ _ exp {ikxAx/2) - exp {-ikxAx/2) _ ^ sin {kxAx/2) 

iAx ~ {kxAx/2) ■ ' 

Our temporal discretization ensures that the discrete velocities are discretely divergence free, that 
is, {viv\) = to within the tolerance of the linear solvers used for the Stokes system. For a perfect 
scheme, the dimensionless structure factor 

and analogously Si, (in three dimensions) would be unity for all wavenumbers, while 5'i ~ (^2^^3) 
would be zero. 

Note that for a system at equilibrium, Vc = 0, the linearized velocity equation and the concen- 



tration equation (18) are uncoupled and thus Sc^v = 0. Observe that the same temporal discretiza- 
tion is used for the velocity equation, projected onto the space of discretely divergence-free vector 
fields consistent with the boundary conditions, and for the concentration equation. Therefore, it is 
sufficient to present here numerical results for only one of the self-correlations s!;^^ and Sc,c- 
In Fig. [2] we show S^"^ as a function of the wavenumber k in three dimensions for a cell Reynolds 
number r = 1 and an advective CFL number a = 0.5 and a = 0.25. Even for the relatively large 
time step, the deviation from unity is less than 5%, and as a — )• it can be shown theoretically 
and observed numerically that the correct covariance is obtained at all wavenumbers. 

Theoretical analysis suggests that the error in the discrete covariance vanishes with time step 
and the background velocity as O(a^) ~ O (V^At^^ for both velocity and concentration [46j. In 
the left panel of Fig. [3] we show the observed relative error in the variance of the discrete velocity as 
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Figure 2: Spectral power of the first solenoidal mode for an incompressible fluid, Sv {kx,ky,kz), as a 
function of the wavenumber (ranging from to tt/Ax along each axes), for a periodic system with 32"^ cells. 
A uniform background flow along the z axis is imposed. The left panel is for a time step a — 0.5, and the 
right for a ~ 0.25. Though not shown, we find that S-i and Sc,c are essentially identical, and both the real 
and imaginary parts of the cross-correlation Si, ' vanish to within statistical accuracy. 

a function of a, confirming the predicted quadratic convergence. As expected, identical results are 
obtained for concentration as well. These numerical results confirm that our spatial discretization 
satisfies discrete fluctuation-dissipation balance and the temporal discretization is weakly second- 
order accurate. 

B. Compressible Solver 

Unlike the incompressible method, which requires complex linear solvers and preconditioners, 
the explicit compressible scheme is very simple and easy to parallelize on Graphics Processing 
Units (GPUs). Our implementation is written in the CUDA programming environment, and is 
three-dimensional with the special case of = 1 cell along the z axes corresponding to a quasi 
two-dimensional system. In our implementation we create one thread per cell, and each thread only 
writes to the memory address associated with its cell and only accesses the memory associated with 
its own and neighboring cells. This avoids concurrent writes and costly synchronizations between 
threads, facilitating efficient execution on the GPU. Further efficiency is gained by using the GPU 
texture unit to perform some of the simple computations such as evaluating the equation of state. 
Our GPU code running in a NVIDIA GeForce GTX 480 is about 4 times faster (using double 
precision) than a compressible CPU-based code [28] running on 32 AMD cores using MPI. Note that 
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Figure 3: (Left) Relative error in the equilibrium variance of velocity (or, equivalently, concentration) for 
several time steps, as obtained using our incompressible code with a background flow velocity Vq = (\/3, 2) /2 
corresponding to cell Reynolds number r =^ %/3/2 in two dimensions, and Vq = (1, 1/3, 1/3) corresponding 
to r = 1 in three dimensions, for a grid of size 32^^ and 32^ cells, respectively. The theoretical order 
of convergence O(Ai^) is shown for comparison. Error bars are on the order of the symbol size. (Right) 
Normalized covariance of the discrete velocities and densities compared to the theoretical expectations, using 
the parameters reported in the caption of Fig. |4] The value reported is the relative error of the variance 
of a variable or the correlation coefficient between pairs of variables, see legend. The theoretical order of 
convergence 0{At^) is shown for comparison. Error bars are indicated but are smaller than the symbol size 
except for the smallest time step. 

with periodic boundary conditions the velocity and the pressure linear systems in the incompressible 
formulation decouple and Fast Fourier Transforms could be used to solve them efficiently. We have 
used this to also implement the incompressible algorithm on a GPU by using the NVIDIA FFT 
library as a Poisson/Helmholtz solver. We emphasize, however, that this approach is applicable 
only to the case of periodic boundary conditions. 

We first examine the equilibrium discrete Fourier spectra of the density and velocity fluctuations 
for a uniform periodic system with an imposed background flow, with similar results observed for 
concentration fiuctuations. In Fig. [4] we show the correlations of density and velocity fluctuations 
as a function of the wavenumber k in three dimensions for a CFL number of a = 0.25. We see that 
self-correlations are close to unity while cross-correlations nearly vanish, as required, with density 
fiuctuations having the largest relative error of 5% for the largest wavenumbers. 

Calculating cross-correlations in real space is complicated by the staggering of the diffent grids. 
We arbitrarily associate one half of the cell faces with the cell center, defining {{Sp) {6vx)) = 
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Figure 4: Normalized static structure factors S'p p (top left), S'„^„^ (top right), Sp^y^ (bottom left) and Sy^^y^ 
(bottom right) for a compressible fluid with physical properties similar to water, for a periodic system with 
30"^ cells. A uniform background flow with velocity vq — (0.2, 0.1, 0.05)ct is imposed and the time step 
corresponds to an acoustic CFL number a = 0.25 and viscous CFL number = 0.017 for shear viscosity 
and = 0.041 for bulk viscosity. 



{6p. 



. Theoretical analysis suggests 



« ^)) and ((foj J 
that the error in the discrete covariance vanishes with time step as O(a^) ~ O (c|iAt^) [?6]. In 
the right panel of Fig. [3] we show the relative error in the discrete covariances as a function of a in 
the presence of a background flow, confirming the predicted cubic convergence. These numerical 
results verify that our spatial discretization satisfies discrete fluctuation-dissipation balance and 
the temporal discretization is weakly third-order accurate. 
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Figure 5: Numerical data (symbols) and theory (lines) for the real part of several dynamic structure factors 
for wavenumber fc — (2, 2, 2) • 27r/L in a cubic periodic box of 30"^ cells and volume L^. Self correlations are 
shown in the left panel, and cross-correlations are shown in the right panel. The imaginary part vanishes 
to within statistical accuracy for the off-diagonal terms. The physical parameters are as reported in the 
caption of Fig. |4] 

1. Dynamic Correlations 



For compressible flow, the dynamics of the fluctuations is affected by the presence of sound 
waves and it is important to verify that the numerical scheme is able to reproduce the temporal 
correlations between the fluctuations of the different pairs of variables. In particular, a good method 
should reproduce the dynamic correlations at small wavenumbers and wave-frequencies correctly 
|28| . Theoretical predictions for the equilibrium covariances of the spatio-temporal specta of the 
fluctuating fields, usually referred to as dynamic structure factors, are easily obtained by solving 
the equations (6|7) in the Fourier wavevector-frequency {k,uj) domain and averaging over the 
fluctuations of the stochastic forcing [17J. The density-density dynamic structure factor Sp^p{k,u}) 
is accessible experimentally via light scattering measurements, and for isothermal flow it exhibits 
two symmetric Brilloin peaks at a; w ztc-rA;. The velocity components exhibit an additional central 
Rayleigh peak at uj = due to the viscous dissipation. As the fluid becomes less compressible (i.e., 
the speed of sound increases) , there is an increasing separation of time-scales between the side and 
central spectral peaks, showing the familiar numerical stiffness of the compressible Navier-Stokes 
equations. 

In Fig. [5] we compare the theoretical to the numerical dynamic structure factors for one of 
the smallest resolved wavenumbers, and observe very good agreement. Note that unlike static 
correlations, dynamic correlations are subject to discretization artifacts for larger wavenumbers. 
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even as At — )■ [28]. Specifically, the positions and widths of the various peaks are set by the 
effective wavevector k rather than the true wavevector k, as given for the standard second-order 



discretization of diffusion in (57). 



V. GIANT FLUCTUATIONS 



As a non-trivial application of our staggered schemes for fluctuating hydrodynamics, we perform 
the first incompressible computer simulations of diffusive mixing in microgravity, recently studied 
experimentally aboard a satellite in orbit around the Earth [12j. The experimental data presented in 
Ref. [T5J shows good agreement with theoretical predictions, however, various over-simplifications 
are made in the theory, notably, only the solenoidal velocity mode with the largest wavelength 
is considered. Numerical simulations allow for a more detailed comparison of experimental data 
with fluctuating hydrodynamics, at least within the applicability of the physical approximations 
discussed in Section [TAl 

The experimental conflguration consists of a dilute solution of polystyrene in toluene, confined 
between two parallel transparent plates that are a distance h = 1mm apart. A steady temperature 
gradient VT = /ST /h is imposed along the y axes via the plates. The weak temperature gradient 
leads to a strong concentration gradient Vc = c5tVT due to the Soret effect, giving rise to an 
exponential steady-state concentration profile c{y). Quantitative shadowgraphy is used to observe 
and measure the strength of the fluctuations in the concentration around c via the change in the 
refraction index. The observed light intensity, once corrected for the optical transfer function of 
the equipment, is proportional to the intensity of the fluctuations in the concentration averaged 
along the gradient, 

c±{x,z) = h ^ c{x,y,z)dy. 

Jy=0 

The main physical parameters we employed in our simulations are summarized in Table |I} Addi- 
tional details of the experimental setup and parameters are given in Ref. [12j. 

The large speed of sound in toluene makes the compressible equations very stiff at the length 
scales of the experimental system. It is usually argued that compressibility does not affect the 
concentration fluctuations [17j . Solving the compressible equations in the presence of a concen- 
tration gradient conflrms that, as long as there is a large separation of time scales between the 
acoustic and diffusive dynamics, the presence of sound waves does not affect the concentration 
fluctuations. In our compressible simulations, we artificially decrease the speed of sound many- 
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Parameter 


Value 


Notes 


P 


0.86 gr/cm^ 


On average only if compressible 


x{v + x) 


1.2 • 10~^ cm^/s^ 


Kept constant in all runs 


V 


Variable Sc — v/x 


Physical value v = 6.07 • 10^'' cm^/s 


X 


Variable Sc = vjx 


Physical value % = 1.97 • 10^^ cm^/s 


C 





None for incompressible 




4.18 • gr cmVs 


Corresponds to T = 303 K 


M 


1.51 • 10^20 gj. 


Not important for results 


St 


0.0649 K"^ 


Enters only via S'tVT 


Co 


0.018 


On average only if nonperiodic 


CT 


1.11 cm/s 


Physical value ct ~ 1.3 • 10^ cm/s 



Table I: Summary of parameters used in the simulations of giant fluctuations in zero gravity. 



fold and set the cell Reynolds number to r = CTjiy^x) > 10. Numerical results show that this 
is sufficient to approach the limit r — )• cxd to within the statistical accuracy of our results. This 
decrease in ct corresponds to making the mass of the toluene molecules much larger than the 
mass of the polystyrene macromolecules themselves, which is of course physically very unrealistic. 
One can think of our compressible simulations of giant fluctuations in microgravity as an artificial 
compressibility method for solving the incompressible equations. 

In the actual experiments reported in Ref. [12j, concentration diffusion is much slower than 
momentum diffusion, corresponding to Schmidt number S'c = z^/x ^ 3 • 10^. This level of stiffness 
makes direct simulation of the temporal dynamics of the fluctuations infeasible, as long averaging is 
needed to obtain accurate steady-state spectra, especially for small wavenumbers. However, as far 



as the nonequilibrium static correlations are concerned, we see from (27) that the crucial quantity 
is + x) = ("5 + rather than x and v individually. Therefore, we can artificially increase 

X and decrease to reduce s, keeping s 1 and (s + l)x^ fixed. In the linearized case, it can be 
proven more formally that there exists a limiting stochastic process for the concentration as s — ?• oo 
so long as sx^ is kept constant (E. Vanden-Eijnden, private communication). In fact, artificially 
decreasing the Schmidt number while keeping sx^ fixed can be seen as an instance of the seamless 
multiscale method presented in Ref. 
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A. Approximate Theory 

For large wavenumbers the influence of the boundaries can be neglected and the periodic theory 



presented in Section |IIB 1| applied. In order to demonstrate the importance of the boundaries, 
and also to test the code by comparing to the periodic theory, we have implemented a model in 
which qualitatively similar giant concentration fluctuations appear even though the macroscopic 
concentration profile is uniform, c(y) = cq. Numerically, this sort of quasi-periodic model is 
implemented by using periodic boundary conditions but adding an additional source term —v ■ Vc 



in the concentration equation, as in (18). This term mimics our skew-adjoint discretization of the 
advection by the fluctuating velocities 



and is conservative when integrated over the whole domain. Note that in this quasi-periodic setup 
Vc is simply an externally-imposed quantity unrelated to the actual mean concentration profile. We 
emphasize that these quasi-periodic simulations are used only for testing and theoretical analysis of 
the problem, and not for comparison with the experimental results. In the simulations with physical 
boundaries and in the experiments the concentration profile is exponential rather than linear. For 
the purposes of constructing a quasi-periodic approximation we take the effective concentration 
gradient to be Vc w Ac/h, where Ac is the difference in concentration near the two boundaries. 
For periodic systems, the spectrum of the fluctuations of c± can be obtained from the full three- 



dimensional spectrum (27) by setting ky = /cy = 0. For the specific parameters in question the 
equilibrium fluctuations in concentration are negligible even at the largest resolved wavenumbers. 
When discretization artifacts are taken into account, the quasi-periodic theoretical prediction for 
the experimentally-observed spectrum becomes 

5c^P (A:., k.) = ( (6c A (6cX) = /"^^ ~, (Vc)^ , (58) 



where k^ = (k^. + k'l\ and tilde denotes the effective wavenumber (57). Imposing no-slip condi- 



tions for the fluctuating velocities makes the theory substantially more complicated. A single-mode 
approximation for the velocities is made in Ref. [i_62J in order to obtain a closed-form expression 
for the spectrum of concentration fluctuations in a non-periodic system S^^. For a small Lewis 
number and without gravity it is found that 

^p(fcl) " ^^^^^^ = gi + 24.6^1 + 500.5' ^^^^ 
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Figure 6: Snapshots of the concentration c± in the plane perpendicular to the gradient Vc, at times O.lro, 
To, and 5to after the gradient is established. The thickness of the sample (perpendicular to the page) is 
one quarter of the lateral extents of the system, h — Ly ~ Lx/4:, and sets the scale of the steady-state 
fluctuations. Compare to the experimental snapshots shown in Fig. 1 of Ref. [12j . 



where q± = hk± is a non-dimensionalized wavenumber. 



The Galerkin function G given by ( 59 ) reflects the physical intuition that the no-slip condition 



suppresses fluctuations at scales larger than the distance between the physical boundaries 
After the concentration gradient is established, "giant" [32] concentration fluctuations evolve with 
a typical time scale of tq = h? /{ir'^x) ~ 1000s, until a steady state is reached in which the typical 
length scale of the concentration fluctuations is set by the finite extent of the domain. This is 
illustrated in Fig. [g] via snapshots of c_\_{x, z; t) taken at several points in time after starting with 
no concentration fluctuations at time t = 0. 



B. Simulations and Results 



In our simulations, the plates are represented by no-slip boundaries at y = and y = h, and 
periodic boundaries are imposed along the x and z axis to mimic the large extents of the system in 
the directions perpendicular to the gradient. A Robin boundary condition is used for concentration 
at the physical boundary, 



dc 
dn 



-c{n ■ Vs) , 



ensuring that the normal component of the concentration flux vanishes at a physical boundary. 
The stochastic concentration flux also vanishes at the boundary as for Dirichlet boundaries, since 
the Soret term does not affect fluctuation-dissipation balance. In the codes the boundary condition 



41 



is imposed by setting the concentration in a ghost ceh to 

2 ± VsAy 

where c„ is the value in the neighboring ceh in the interior of the computational domain, and the 
sign depends on whether the ghost cell is at the low or high end of the y axes. The boundary 
condition is imposed explicitly, which leads to non-conservation of the total concentration when 
a semi-implicit method is used for the diffusive terms in the concentration equation. This can 
be corrected by implementing the boundary condition implicitly or using an explicit method for 
concentration; however, we do not do either since the observed change in the average concentration 
is small for the specific parameters we use. 

Using the incompressible formulation allows for a much larger time step, not only because of 
the lack of acoustics, but also because of the implicit temporal discretization of the viscous terms 
in the momentum equations. However, it is important to remember that a time step of our GPU- 
parallelized compressible code takes much less computing than a time step of the incompressible 
code. Nevertheless, we are able to study larger system sizes in three dimensions using the in- 



compressible algorithm. In the incompressible sijmulations, we used (36) for the velocity equation 
in order to avoid unnecessary projections. Because of the explicit handling of the concentration 
boundary conditions, we employed a predictor-corrector algorithm for the concentration equation. 



in which both the predictor and the corrector stages have the form (37). 

In Fig. [7] we show numerical results for the steady-state spectrum of the discrete concentration 
field averaged along the y-axes, in two (left panel) and in three dimensions (right panel), for both 
bulk (quasi-periodic) and finite (non-periodic) systems. In order to compare with the theoretical 



predictions (58) and (59) most directly, we plot the ratio of the observed to the predicted spectrum. 
This choice of normalization not only emphasizes any mismatch with the theory, but also eliminates 
the power-law {kj_'^) divergence and makes it easier to average over nearby wavenumbers k± and 
also estimate error bars^. For the runs reported in Fig. [7] we applied the largest concentration 
(temperature) gradient (AT = 17. 4K) used in the experiments [12]; we have verified that the 
non-equilibrium concentration fiuctuations scale as the square of the gradient. 

Both panels in Fig. [Tjshow an excellent agreement between the theory (58) and the numerical 
results for quasi-periodic systems. This shows that correcting for the spatial discretization artifacts 
by replacing k± with k± accounts for most of the discretization error. For the compressible runs, we 



^ Note, however, that the most rehable error bars are obtained by averaging over many uncorrelated runs started 
with difTerent random number seeds. 
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Figure 7: Ratio between the numerical and theoretical discrete spectrum of concentration projected along 
the y axes, as a function of the normalized wavenumber q± = k±h. For all runs Ny = 32 cubic hydrodynamic 
cells along the y axes were used, and all systems have aspect ratio N^/Ny — N^/Ny = 4. Error bars are 
indicated for some of the curves to indicate the typical level of statistical uncertainty. (Left) Two dimensions, 
for both compressible and incompressible fluids (see legend), with either periodic boundary conditions (empty 
symbols) or physical boundaries (solid symbols) imposed at the y-boundaries, for several Schmidt numbers 
s — vjx- (Right) Three dimensions, with same symbols as left panel), along with arbitrarily normalized 
experimental data |T2] (see legend) corresponding to s w 3 • lO'^ (courtesy of A. Vailati). 



use a relatively small time step, a = 0.2, leading to temporal discretization errors that are smaller 
than the statistical accuracy except at the largest wavenumbers. Our semi-implicit discretization 
of the incompressible equations gives the correct static covariance of the concentration for all time 
step sizes. Based on the analysis presented in Appendix |XJ the majority of the incompressible 
simulations employ a time step corresponding to a viscous CFL number /? = 1 or /3 = 2, with a 
few of the largest systems run at /3 = 5 to resolve the smaller wavenumbers better. 

In the left panel of Fig. [7] we compare results from two-dimensional compressible and in- 
compressible simulations and find excellent agreement. For non-periodic systems the single- mode 



Galerkin theory (59 ) is not exact and the theory visibly over-predicts the concentration fluctuations 



for smaller wavenumbers in both two and three dimensions. We observe only a partial overlap of 
the data for different Schmidt numbers s = vjx for smaller wavenumbers, although the difference 
between s = 10 and s = 20 is relatively small. 

In three dimensions, we rely on the incompressible code in order to reach time scales necessary 
to obtain sufficiently accurate steady-state averages for large Schmidt numbers. In the right panel 
of Fig. [7] we compare numerical results for quasi-periodic and non-periodic compressible and 
incompressible systems to the theoretical predictions and also to experimental data from Ref. 



43 



|12j (A. Vailati, private communication). While the numerical data do not match the experiments 
precisely at the smallest wavenumbers, a more careful comparison is at present not possible. Firstly, 
the boundary conditions affect the small wavenumbers strongly, and our use of periodic boundary 
conditions in x and z directions does not match the experimental setup. The experimental data 
has substantial measurement uncertainties, and is presently normalized by an arbitrary pre-factor. 
Within this arbitrary normalization, our numerical results seem to be in good agreement with 
the experimental observations over the whole range of experimentally-accessible wavenumbers, 
and the agreement at small wavenumbers improves as the Schmidt number of the simulations 
increases. The actual magnitude of the macroscopic non-equilibrium fluctuations in c_l is given 
by the integral of the structure factor S^^ over all wavenumbers k^. Numerically we observe 
fluctuations ^((5c_l)^^ /(?[_ ~ 3 • 10^^, which is consistent with experimental estimates (A. Vailati, 
private communication) . 



VI. CONCLUSIONS 



We have presented spatio-temporal discretizations of the equations of fluctuating hydrodynam- 
ics for both compressible and incompressible mixtures of dynamically-identical isothermal fluids. 
As proposed by some of us in Ref. |28], we judge the weak accuracy of the schemes by their ability 
to reproduce the equilibrium covariances of the fluctuating variables. In particular, for small time 
steps the spatial discretization ensures that each mode is equally forced and dissipated in agree- 
ment with the fluctuation-dissipation balance principle satisfied by the continuum equations. A 
crucial ingredient of this discrete fluctuation-dissipation balance is the use of a discrete Laplacian 
L = —DD* for the dissipative fluxes, where D is a conservative discrete divergence, with a suitable 
correction to both the Laplacian stencil and the stochastic fluxes at physical boundaries. Further- 
more, we utilize a centered skew-adjoint discretization of advection which does not additionally 
dissipate or force the fluctuations, as previously employed in long-time simulations of turbulent 
flow, where it is also crucial to ensure conservation and avoid artificial dissipation |77j . 

For the compressible equations, our spatio-temporal discretization is closely based on the col- 
located scheme proposed by some of us in Ref. [28], except that here we employ a staggered 
velocity grid. It is important to point that out the difference between a collocated scheme, in which 
the fluid variables are cell-centered but the stochastic fluxes are face-centered (staggered), as de- 
scribed in Ref. |28j, and a centered scheme where all quantities are cell-centered. Several authors 
p6l [27] have already noted that centered schemes lead to a Laplacian that decouples neighboring 
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cells, which is problematic in the context of fluctuating hydrodynamics. We emphasize however 
that these problems are not shared by collocated schemes for compressible fluids, for which the 
Laplacian L = —DD* has the usual compact 2d+l stencil, where d is the dimensionality |2H]- Dis- 
cretizations in which all conserved quantities are collocated may be preferred over staggered ones in 
particle-continuum hybrids or more generally, in conservative discretizations for non-uniform 
grids. 

A staggered grid arrangement, however, has a distinct advantage for incompressible flow. 
Namely, the use of a staggered grid simplifies the construction of a robust idempotent discrete 
projection P = / -|- D*L^^D that maintains discrete fluctuation-dissipation at all wavenumbers. 
In the temporal discretization employed here, based on prior work by one of us p34|, this projection 
is used as a preconditioner for solving the Stokes equations for the pressure and velocities at the 
next time step. For periodic systems the method becomes equivalent to a classical Crank-Nicolson- 
based projection method, while at the same time avoiding the appearance of artificial pressure 
modes in the presence of physical boundaries [72] . 

The numerical results presented in Section |V] verify that our numerical simulations model ex- 
perimental measurements of giant fluctuations [12] during diffusive mixing of fluids faithfully. The 
numerical simulations give access to a lot more data than experimentally measurable. For example, 
the spectrum of concentration fluctuations in the x — z plane can be computed for planes (slices) as 
the distance from the boundaries is varied, giving a more complete picture of the three dimensional 
spatial correlations of the nonequilibrium fluctuations. We defer a more detailed analysis, including 
a study of temporal correlations, to future work. 

The compressible solver we developed utilizes modern GPUs for accelerating the computations. 
In the future we will investigate the use of GPUs for the incompressible equations, starting with 
simple FFT-based solvers for periodic systems. For grid sizes that are much larger than molecular 
scales, the stability restriction of explicit compressible solvers becomes severe and it becomes nec- 
essary to eliminate sound waves from the equations by employing the low Mach number limit. A 
challenge that remains to be addressed in future work is the design of zero Mach number methods 
|48j for solving the variable-density equations of fluctuating hydrodynamics, as necessary when 
modeling mixtures of miscible fluids with different densities. This would enable computational 
modeling of the effects of buoyancy (gravity) in experimental studies of the giant fluctuation phe- 
nomenon performed on Earth |39t 1^ . 
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Appendix A: Implicit Midpoint Rule as a Gibbs Sampler 

We consider here numerical methods for the general additive-noise linear SDE 

— =Ax + KW {£) , (Al) 

where W(t) denotes white noise. If the eigenvalues of A have negative real parts, the long-time 
dynamics tends to a Gaussian equilibrium distribution 

Peq(^) = ^-^expf-^^^V (A2) 



2 

where the covariance matrix S is the solution to the linear system [see, for example, Eq. (30) in 
[28] or Eq. (3.10) in [6l]] 

AS + SA* = -KK\ (A3) 
If one is only interested in calculating steady-state observables (expectation values), then a numer- 



ical method for solving (Al) needs to only sample the equilibrium Gibbs distribution (A2), without 



having to approximate the correct dynamics. 



The implicit midpoint rule or Crank-Nicolson discretization that we employed in Section III A 2 

/t" -I- T^+l \ 

= + ^ ^ At^/^KW, (A4) 
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can be seen as a Markov Chain Monte Carlo (MCMC) algorithm for sampling from the distribution 



(A2). This sampling is exact, that is the equilibrium distribution of the chain (A4) is exactly (A2). 



This important fact can be shown using the techniques described in Ref. [S^, but here we present 
an alternative derivation. 

A well-known MCMC algorithm for sampling the Gibbs distribution is the Metropolis-Hastings 
algorithm. In this algorithm, one treats x^'^^ as a trial or proposal move that is then to be accepted 
with probability 



Peg Prev (3^"+^ ^ X^) 



where Pforw is the transition probability for the chain (A4) and Prev is the transition probability 



for the time-reversed chain (this important distinction ensures strict time reversibility of the chain 
with respect to the equilibrium distribution). Explicitly, 



Prev 



P 



forw 



X 



C exp 
C exp 



W 



W 



where the reverse step noise W is the solution to the equation (here the adjoint of A appears 
because of time reversal) 



+ A" 



X + X 



n+l 



A tedious but straightforward matrix calculation shows that the acceptance probability a = 1, 
that is, no rejection is necessary for the implicit midpoint rule to sample the correct equilibrium 
distribution, regardless of the time step At. 

The calculation of a is simple to do if a Fourier transform is used to diagonalize the hydrody- 



namic equations [see Eq. (19)] to obtain a system of scalar SDEs with complex coefficients. For the 



stochastic advection-diffusion equation (38) with v = vq, which is a good model for more general 
hydrodynamic equations, 



A 



-a + bi, K = K 



2a, and S = S = 1, 



(A5) 



with a = and h = —kvo, where k is the wavenumber. 

While the time step At can be chosen arbitrarily without biasing the sampling, the optimal 
choice is the one that minimizes the variance of the Monte Carlo estimate of the observable of 
interest. In the simulations of giant fluctuation experiments, the observable of interest is the 
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covariance (spectrum) of the fluctuations S = {xx*) . The variance of the Monte Carlo estimate of 
S is proportional to the autocorrelation time r of = a;"" (a;")*, which is itself proportional to 
the sum of the autocorrelation function of S"' [5S]. Focusing on the scalar SODE A5, we get the 
autocorrelation time 



E 

n=0 



n=0 



1 1 aAt bAt 



For the purely diffusive equation, vq = 0, the statistical accuracy for a fixed number of time steps 
is proportional to 



~ 4 + 4A:2/3 + ^4^2' 

where /3 = vAtjAx^ is the viscous CFL number and k = kAx is the dimensionless wavenumber. 
Note that r^^ ~ (3k'^ for small fc, so increasing the time step improves the sampling. However, for 
large k increasing the time step reduces the statistical accuracy (this is related to the fact that the 
Crank-Nicolson algorithm is A-stable but it is not L-stable), . The wavenumber 

with highest statistical accuracy fcopt depends on the time step, /3feopt — 2, or, alternatively, the 
optimal choice of time step depends on the wavenumber of most interest. For the type of problems 
we studied in this work the spectrum of the fluctuations has power-law tails ~ k~^ and therefore 
all wavenumbers are important. Using /? ~ 2 produces a good coverage of all of the wavenumbers. 



Appendix B: Fluctuation-Dissipation Balance for Incompressible Flow 

Discrete fluctuation-dissipation balance is affected by the presence of an incompressibility con- 
straint. The spatially discretized velocity equation linearized around a stationary equilibrium state 
has the form, omitting unimportant constants in the noise amplitude, 



dtv 



vLv + V2uDWy 



(Bl) 



where we used a non-symmetric stochastic stress tensor since the symmetry does not affect the 
results presented here. The steady-state covariance of the velocities = (vv*) is determined from 



the fluctuation-dissipation balance condition (A3) with A = uFL and K = v2i^PZ), giving 



FLS^ + S^L-'F* = -IFDD-'F*. (B2) 
The fluctuation-dissipation balance condition for the simple advection-diffusion equation, 

L + L* = DG + {DGY = -2DD*, 
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as stated in (45) with all of 



implies that = P is the solution to (B2) if P is self-adjoint, P* = 
the constants included. 

The above analysis was does not account for the temporal discretization. For small timesteps, 



our temporal discretization of (Bl) behaves like a projected Euler-Maruyama method, 

' [i?" + uLvAt + V2ijAtDWv . 



An important difference with the continuum equation (Bl) is that the velocity in the previous 



time step is also projected, i.e., the increment of O (At) is added to Pd" and not to d". If P is 



idempotent. 



just as the continuum projection operator is, then subsequent applications of 



the projection operator do not matter since v"^ is already discretely divergence free, = v"'. In 
the literature on projection methods idempotent projections are called exact projections. 



The above considerations lead to the conclusion that S„ 



'if] 



' and ^ 



Both of these 



conditions are met by the MAC discrete projection operator P = / — D* (DD*)^^ D, which shows 
that our spatio-temporal discretization gives velocity fluctuations that have the correct covariance 
( 45 ) . A straightforward extension of the analysis in Appendix [A] shows that the Crank-Nicolson 



temporal discretization (36) gives the correct equilibrium velocity covariance for any timestep size. 



not just for small time steps. Further details will be presented in future publications [l6] 
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